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1. Introduction 

1.1. Background 

The contemporary mesoscopic theoretical description of multiple scattering of waves 
in random media is most often based on the linear radiative transport equation 
(RTE) [1]. Unfortunately, the RTE is notoriously difficult to solve, even in the case 
of constant absorption and scattering coefficients. The known analytical solutions are 
few and of little practical importance. Yet, there is a growing need for accurate and 
computationally efficient solutions to the RTE in many fields of applied and fundamental 
science. For example, in optical tomography of biological tissues [2,3], the use of the 
RTE is frequently required to accurately describe propagation of multiply scattered 
light. This is especially true in close proximity to sources or boundaries [4], or in 
regions with high absorption and low scattering [5,6]. Accordingly, significant effort has 
been devoted to developing and refining efficient approximate and numerical methods 
for solving the RTE. In particular, recently explored approaches have been based on 
the discrete ordinate method [7-9], cumulant expansion [10,11], modifications of the 
Ambarzumian's method [12,13], and different levels of the Pl approximation [14,15]. 
Algorithms for inversion of the RTE have also been proposed [16-19]. 

The discrete ordinate method (see [20] for a detailed description) is, perhaps, the 
most common approach due to its simplicity and generality. An alternative to the 
discrete ordinate method is the method of spherical harmonics, often referred to as the 
Pl approximation, in cases with special symmetry. This approach has the advantage 
of expressing the angular dependence of the specific intensity in a basis of analytical 
functions | rather than in the completely local basis of discrete ordinates. In particular, 
in the case of cylindrical symmetry (one-dimensional propagation), a very effective 
solution based on a contuned fraction expansion can be obtained [21]. However, when no 
special symmetry is present in the problem, the method of spherical harmonics can be 
carried out in practice only to very low orders. In a recent paper [22] , we have suggested 
a modification of the standard method of spherical harmonics. The modification is based 
on expanding the angular part of each Fourier component of the specific intensity in 
the basis of spherical functions defined in a reference frame whose ^-axis is aligned 
with the direction of the Fourier wave vector. This approach resulted in significant 
mathematical simplifications and was referred to as the modified method of spherical 
harmonics in [22]. Here we find it more appropriate to call it the method of rotated 
reference frames (MRRF). 

In Ref. [22], the derivation of the RTE Green's function by the MRRF was only 
briefly sketched in infinite-medium and numerical examples were limited to a few simple 
cases with spherical symmetry. Here we give the full mathematical details of the 
derivation and discuss the mathematical properties of the solutions obtained, derive 

I Here the term "analytical" means "expressed in terms of well-characterized functions through explicit 
formulas" , not necessarily analytic in the Cauchy-Riemann sense. 
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plane-wave decomposition of the Green's function, and generalize the MRRF to the 
case of planar boundaries. We also provide extensive numerical examples for cases with 
no special symmetry. The paper is organized as follows. In Section II. 2| we introduce 
the RTE and basic notations, and explain why the use of rotated reference frames is 
beneficial. In Section El we define spherical functions in rotated reference frames. In 
Section El we apply the MRRF to the derivation of the Green's function. In particular, 
the Green's function in the Fourier representation is given in Section ITTl Mathematical 
properties of the solutions are discussed in Section 021 Different representations for the 
Green's function in real space are given in Section IH7S1 A plane-wave decomposition of 
the Green's function is derived in Section 03 In Section ITHl we introduce evanescent 
modes of the homogeneous RTE. These modes are important mathematical constructs 
which can be used for solving the RTE in the presence of planar boundaries, as is 
shown in Section IH.fil Section |31 contains numerical examples of applying the MRRF 
to calculating the Green's function in infinite space. Finally, Section El contains a 
discussion. 



1.2. RTE and the conventional method of spherical harmonics 

The RTE describes the propagation of the specific intensity /(r, s), at the spatial point 
r and flowing in the direction specified by the unit vector s, in a medium characterized 
by absorption and scattering coefficients fia and Hs, and has the form 

s ■ V/ + (/ia + = fisAI + e . (1) 

Here e = e{r, s) is the source and A is the scattering operator defined by 

AI{r,s) = J A{s,s')I{r,s')d^s' . (2) 

The phase function A{s, s') is normalized according to the condition J A{s, s')d'^s' = 1. 
We also assume that it depends only on the angle between s and s': A{s, s') = /(s ■ s'). 
This fundamental assumption is often used and corresponds to scattering by spherically 
symmetric particles. 

In the conventional method of spherical harmonics, all angle-dependent quantities 
are expanded in the basis of spherical harmonics defined in the laboratory frame [23]: 



/(r,s) =5^/,„(r)rz„(s) , (3) 

bn 

e{r,s) = ^eim{r)Yirr,{s) , (4) 

Im 

A{s,s') = j2^iyu^)yLi^') ■ (5) 



Im 



In particular, truncating the above series at / = 1 leads to the well-known diffusion 
approximation to the RTE [23]. In the more general case, substituting expansions 
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(0) into the RTE ([T]), multiplying the resulting equations by Yjf,^,(s) and integrating 
over s leads to the following system of equations for //^(r): 



where R^°'^ = J SaYi^(s)Yii.in'{s)d'^s {a = x,y,z) are matrices whose explicit form is 
given in Ref. [23] and 



This system of partial differential equations must be solved for / = 0, 1, 2, ... , /max and 
m = — /, ...,/, where /max is the truncation order of the expansion 

In a classic text, Case and Zweifel wrote concerning the system of equations (jHl): 
"This rather awe-inspiring set of equations . . . has perhaps only academic interest" 
(Ref. [23], p. 219). We note that the root of the difficulty is not that the matrices 

are dense (in fact, they only couple coefficients with m' = m, /' = / ± 1 for a = z 
and m' = mdzl, /' = /±1 for a = x,y) or non-commuting (in fact, it is easy to verify that 
all R^"^ commute). The difficulty is that these matrices operate on the spatial derivatives 
of Iim taken along different directions. Thus, by viewing the set of three matrices R^"^ 
as a three-dimensional vector of matrices R, and using the Fourier representation for 
Iim, we can rewrite the term in the square bracket of (jH)) as ik • Him^i'm'h'm'O^)- It can 
be seen that the matrix k ■ R depends explicitly on the direction and length of k. (See 
similar formulation in Ref. [24].) 

The method of rotated reference frames (MRRF), similar to the conventional 
method of spherical harmonics, does not lead to the separation of spatial and angular 
variables which is impossible for the RTE. However, by choosing a different k-dependent 
angular basis, we replace the dot product k ■ R by an expression of the type kR, where 
k = |k| is a scalar and i? is a single k-independent block-diagonal matrix. It is shown 
below that, given generalized eigenvectors and eigenvalues of R which must be computed 
numerically, the solution can be obtained in terms of analytical functions of spatial and 
angular variables. § 

2. Spherical functions in rotated frames 

The ordinary spherical harmonics Yim{0, are functions of two polar angles in a fixed 
(laboratory) reference frame. Equivalently, we can view them as functions of a unit 
vector, s. In this case, 6 and ip are the polar angles of s in the laboratory frame. 
More generally, both the orientation of the reference frame and the direction of s can 

§ In principle, it should be also possible to use the fact that all matrices i?^"'' in ^ commute and, 
hence, have the same set of eigenvectors, to solve by diagonalizing just one k-independent matrix 
and analytically inverting the Fourier transform, and thus avoid the use of rotated reference frames. 
This approach has some advantages and difhculties associated with it, and to the best of our knowledge, 
has not been explored so far. If successful, it should lead to the same solutions as described below. 




+ O-Jlm — 



(6) 



(7) 
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Figure 1. Illustration of the rotated reference frame 



vary. We will need to define spherical functions of a unit vector s in a reference frame 
whose z-axis coincides with the direction of a given unit vector k. Obviously, there are 
infinitely many such reference frames. To define one uniquely, it is sufficient to consider 
a rotation of the laboratory frame with the following three Euler angles: a = ip^, P = 6^ 
and 7 = 0, where 9^ and ip^ are the polar angles of k in the laboratory frame. The 
transformation from the laboratory frame (x, y, z) to the rotated frame (x', y' , z') is 
illustrated in Fig. ^ We denote spherical functions of s in the reference frame defined 
by the above transformation by k). They can be expressed as linear combinations 

of the spherical functions defined in the original (laboratory) frame according to 

I 

YU^- k) = YU^- z') = Dl,^{ip^, Q)Yi„S: z) , (8) 

m'=-l 

where 

DLm' (a> 7) = exp(-zma)dj„„, exp(-im'7) (9) 

are the Wigner D-functions; the explicit form of d^^^iiP) is given, for example, in 
Ref. [25]. 
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It is important to note that the expansion of the scattering kernel into the spherical 
functions Yjm(s; k) is independent of the direction of k: 

A{s,s')=J2AiYU^;k)YC^{s';k). (10) 

Im 

Here the expansion coefficients Ai are independent of k and are the same as in (jS)). This 
fact follows from the rotational invariance of the scalar product. 

3. Theory 

3.1. Green's function in the Fourier representation 

By definition, the Green's function G{r, s; fq, Sq) satisfies RTE (PJ with the delta-source 
e = 6{r — ro)5(s — Sq). We will refer to fq, Sq and r, s as the location and direction of 
the source and detector, respectively. In infinite isotropic space, the Green's function 
can be written in the following general form: 



G{r, s; ro, Sq) - 
d^k 



^ J (27r)3 

lm,l'm' ^ ' 



exp[«k-(r-ro)]l^™(s;k)(/m|ir(A;)|/W)F;^,(so;k) . (11) 



Here Kifz) is an unknown operator. The reciprocity of the Green's function, 
G(r, s; ro, So) = G'(ro, — So; r, — s), together with the fact that G is real, imply the 
following symmetry property of K: (l'm'\K\lm) = (— l)'^''(/m|if|/'m')*. This can be 
also written as VK'^V = K, where V is the coordinate inversion operator with matrix 
elements {lm\V\l'm') = {—lYdu'Smm' and f denotes Hermitian conjugation. Thus, it can 
be seen that K is not a Hermitian operator. 

Substituting (|TT|l into ((H) and using the orthogonality properties of the spherical 
functions, we arrive at the following operator equation for K{k): 

{ikR + ^)K{k) = 1 . (12) 

The matrices R and S are defined by 

{lm\R\l'm') = y"(s ■ ^)Yi*mi^; ^)Yi'n,is; k)d''s 

= 6mm' [blmSl'=l-l + &«+l,m^«'=«+l] , (13) 

km = v/(/2-m2)/(4P-l) , (14) 

(/m|S|/W) = Cri6u>6mm' , (15) 

where cx; is given by ((Tj). The formal solution to (fT^ can be written as 

K{k) = S{l+ikW)-^S , (16) 
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where S = 1/vS and W = SRS. Note that {lm\S\l'm') = ?iiv^mm,' I \f^i exists because 
Gi > 0, which follows from the inequalities Ai < 1 and /ia > ||. Similarly to R, 
W is a. real symmetric matrix. Therefore, we can use the spectral theorem to express 
{l + ikW)~^ in terms of the eigenvectors and eigenvalues of W, \ipii) and A^, respectively. 
This immediately leads to the following expression for K{k): 

Given the set of eigenvectors and eigenvalues, which can be found by numerical 
diagonalization of W , the above formula solves the problem in Fourier space. Since 
the components of in the \lm) basis are purely real, it can be seen that K is 
symmetric. 



3.2. Mathematical properties of the solution 

3.2.1. Block structure ofW 

First, we note that W is block-diagonal: {lm\W\l'm') = 5mm'(^|-B(m)|/'). Below, we 
will label different blocks B{M) (M = 0,±1,±2,...) by the capital letter M. The 
matrix elements of B{M) are given by 



{l\B{M)\l') = l3i{M)5v=i.i+Pi+i{M)5v=i+i , > |M| , (18) 

A(M) = hiM/^/^m~i ■ (19) 

Obviously, to find all eigenvalues and eigenvectors of W , it is sufficient to diagionalize 
each block separately. This task is further simplified because all blocks B{M) are 
tridiagonal. We denote eigenvectors of a block B{M) by |0„(M)). Then the eigenvector 
of the full matrix W with the same eigenvalue is obtained according to 

{lm\ijMn) = SmM{l\MM)) . (20) 

The corresponding eigenvalue is denoted by Xmti, where we have introduced a composite 
index (M,n). Note that B{M) = B{-M). 



3.2.2. Symmetry properties of the eigenvectors 

The property VK'^V = K and Eq. ^ imply that VWV = -W. Thus, W is odd 
with respect to coordinate inversion. In particular, if lip) is an eigenvector of W with 
the eigenvalue A, then lip) = V\tp) is also an eigenvector of W but with an eigenvalue 
of the opposite sign, A = —A. The complete set of eigenvectors {IV'm) • A* ^ 
where Q is the set of all values of the index fi, can then be equivalently rewritten 
as {\ip^), \ip^) ■ /i G fi"^}, where O"*" is the set of indices that correspond to positive 
eigenvalues A^. The set of indices that correspond to negative eigenvalues can be denoted 
as Q-] then Q = Q+UQ~ and fi+ n fi^ = {0}. 

II Purely scattering media with = can be considered separately. 
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Using these properties, one can transform the summation over all values of /i in 
()17|) to summation over /i G (such sums will be denoted as Yl'^ below). This 
fact facilitates the inverse Fourier transformation (see [Appendix A| ) and solution of the 
boundary- value problem discussed in Section ESI 

3.2.3. Continuous and discrete spectra 

Third, the eigenvalues can belong either to the discrete or continuous spectrum. 
It is easy to see that the spectrum is continuous for |A| < where = fJ'a + fJ^s, and 
discrete for |A| > l/^t- Indeed, consider the three-point recurrence relation that follows 
from the equation W\ip) = Mi')'- 

Pi{m){l-l,m\i;)+pi+i{m){l + l,m\^) = \{lm\^) , I > \m\ . (21) 

In general, it has two types of solutions: polynomial and exponential. Consider the 
asymptotic properties of these solutions. In the limit / — >^ cx) we have Ai ^ 0, ai ^ fit, 
him 1/2 and (3i{m) 1/2/it. The recurrence relation then becomes: 

(/ - 1, m\il)) + (/ + !, m\i)) = 2fit\{lm\ip) . (22) 

The polynomial solutions have the asymptotic form {Imlip) = p'^{\^t)-, where p^{x) 
are general orthogonal polynomials of degree / (not to be confused with the associated 
Legendre functions which solve the recurrence in the particular case /i^ = 0). In 
order for this solution to be an eigenvector of W , it must be bounded. Obviously, this 
requirement is equivalent to |A/i(| < 1. Thus, for every A G [— 1//^*, there is a 

polynomial solution to the three-term recurrence relation that is an eigenvector of W . 

For A outside of the interval [— l//it], polynomial solutions are unbounded 
and, therefore, can not be eigenvectors of W . We then consider exponential solutions 
which behave asymptotically as {Imlip) = (±1)' exp(— p/) where p satisfies the equation 
cosh(p) = ±fit\. In order for this solution to be an eigenvector of W, p must be 
positive. But the above equation has positive roots only when |A| > l/fif Note 
that the exponentially decaying eigenvectors have a finite norm, and, hence, belong 
to the discrete spectrum. Further bounds on the discrete spectrum can be inferred 
from the Gershgorin theorem, which states that, for a fixed M, {XmuI < = 
max;>|M|[/3;(M) -|- Pi^i{M)]. It can be easily verified that tq = A/y/SfXa and tm = l//ia 
for \M\ > 0. 

In numerical computations, the infinite-dimensional matrix W must be truncated 
and the continuous spectrum of W approximated by a discrete spectrum. In this paper 
we treat all eigenvalues as discrete. Thus, for example, the expression (fTTj) contains only 
a sum over discrete modes, although, theoretically, summation over the continuous part 
of the spectrum must be expressed as an integral. Note that an expression involving 
only discrete spectra avails itself more readily to numerical implementation. 
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3.3. Green's function in real space 

The dependence of solution (|T7jl on k is analytical. This allows us to obtain the Green's 
function the coordinate representation by Fourier transform. We substitute ()17|1 into 
the ansatz (fTTj) and express the spherical functions Yim{s; k) and y5/m'(so; k) in terms of 
spherical functions defined in the laboratory frame whose z-axis direction is given by a 
unit vector z according to (jS)),©- The direction of the x- and y-axes of the laboratory 
frame is arbitrary. This leads to the following expression: 

^(r, s; ro, So) = ^ ^ Yi™(s; z)(/m|x(r - fq; z)|ZW)r;„,(so; z) , (23) 

Im I'm' 

where 

X(r; z) = J exp(zk ■ r)V{k; z)Kik)V\k; z) (24) 

Here I'(k; z) = exp{—i(f^^Jz) exp{—i6^^Jy) is the rotation operator whose matrix elements 
are given by the Wigner D-functions, {lm\V(k; z)\l'm') = SiifDl^^,{(p^, 9^^, 0), ip^ and 9^ 
are the polar angles of k in the laboratory frame, and J is the angular momentum 
operator (with h = 1). We note that operators V are unitary and, hence, normal: 

= V'^ . However, T) does not commute with K. The fundamental simplification 
obtained by the MRRF is that V is known analytically while K has a simple form given 
by (fT7|) . In particular, given numerical values of \ip^) and A^, the dependence of K{k) 
on k is also known analytically. 

Below, we consider two different cases. In the first case, the direction of the 
laboratory frame 2:-axis coincides with the direction from the source to the detector, 
namely, z = (r — ro)/|r — ro|. This choice of the angular basis is convenient when the 
source and the detector are always placed on the same line, irrespective of the directions 
s and Sq. In the second case, we choose z = Sq. This approach is useful when the source 
is scanned, e.g., over a two-dimensional plane, but its direction Sq is fixed. This situation 
is typical for optical tomography in the slab geometry. The integral (j24j) for the two 
cases is evaluated in [Appendix A| The result is, in the first case: 



(im|x(r;f)|;W> = •£ (-1)" -tcf^^Cf;:)]^^ 

^ ' M=-l j=0 

^' {lM\i;,){i;,\l'M) ^^ f R\ 

X E ^,.-.1... [j-J , (25) 

Here I = mm{l,l'), kn{x) = —i^hn\ix) is the modified spherical Bessel function of 
the first kind (defined without a factor of 7r/2), Cf^^l^^^^ are the Clebsch-Gordan 
coefficients and denotes summation over only such indices /i that correspond to 
positive eigenvalues A^. It can be seen that x(r; r) is diagonal in m and m', which 
corresponds to the invariance of the Green's function with respect to a simultaneous 
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rotation of the vectors s and Sq around the hne connecting the source and the detector. 
Eq. ()25|) can be further simphfied by expressing the eigenvectors {tp^) in terms of the 
eigenvectors \(f)n{fn)) of the smaller blocks B{m) as discussed in Section 13.2.11 This 
result, as well as a number of special cases, were given in Ref. [22] and are not repeated 
here. 

In the case z = Sq, expression ()23|) contains only matrix elements of x(r; Sq) 
with m' = 0. This follows from the fact that Y5/m/(so;so) = Sm'o 

^{21' + l)/47r. The 

corresponding expression for the matrix elements of Sq) is 



I I 



(/m|x(r;So)|/'0) 



7r(2/ + l)aiau 



X 



fil'M ^\l-l'\+2j,mSr^' 
^l,M,\l-V\+2j,Q^l,m,l'fl 2-^ 



{lM\^^){ij,\l'M) 



\l-l'\+2j 



R 



(26) 



Derivation of the above result is analogous to that for x(r; r); see Appendix A for details. 



3.4- Plane-wave decomposition of the Green's function 

Having in mind further applications of the MRRF to solving boundary value problems, 
we derive the plane-wave decomposition of the Green's function. The latter is defined 
by the two-dimensional Fourier integral 



G(r,s;ro,So) = J] 5Z / 7^ ^^Pf^^ " (/^ " Po)] 

Im I'm' ^ ' 

X Y;„,(s;z)(/m|/s:(q;2;-2;o)|rm')F;„,(so;z) . (27) 

Here z is a selected direction in space which coincides with the 2;-axis of the laboratory 
frame, p is a two-dimensional vector in the x — y plane (r = p + 2;z and p • z = 0) and 
the direction of x- and y-axes is arbitrary. By comparing the above expression to (j21I), 
we find that 

/'^ dk / \ 

exp{ik,z)V{q + zk,; z)K [^/^^TMj V\q + zk,- z) . (28) 

Here P(q -|- zkz] z) should be understood as a function of the polar angles of the vector 
k = q -|- zkz in the laboratory frame. The latter are defined by 

cosO = kz/^q^ + kl , sin^ = q/^q^ + kl . (29) 

Integral (j^Hj) can be evaluated analytically. The following expression for the matrix 
elements of ^(q; z) is derived in [Appendix B| 
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(/m|.(q; z)\l'm') = ^"^^^^"^ " "^^^^^^ [sgn(.)]'+^'+-+-' T T T' 

^ m\=—lm2=—V yu 

x4™J^^(?A,)](/mi|^,)^^^Pl^g!M^ (30) 



where 



QM) = V^' + , (31) 

the complex angles iT{x) are defined by the relations 

cos[ir(x)] = vT+x^ , sin[ir(x)] = —ix , (32) 

and the angle is the polar angle of the two-dimensional vector q in the x — y 
plane. The Wigner d-functions d^^^,{iT) in the above expression are algebraic functions 
of cos(ir) and sin(zr) (an explicit expression in terms of Jacobi polynomials is given 
in [Appendix B| ) . An expression for K(q; 2) in terms of the block eigenvectors |0„(M)) 
which were introduced in Section r3.2.1l is also given in Appendix B Here we note some 
of the symmetry properties of the matrices fi;(q; z). Under inversion of the 2;-axis, we 
have ft(q, —z) = V^n^q, z)Vz, or, in component form, 

{lm\K{q; -z)\l'm') = (-l)'+''+'"+"'(/m|K(q; z)\l'm') . (33) 

Simultaneous inversion of the x- and y-axes (or, equivalently, rotation around the z-axis 
by the angle vr) is expressed as k(— q, z) = Vxyi^{ci, z)Vxy, or, in component form, 

(/m|/t(-q; z)\l'm') = (-l)"^+'"'(/m|K(q; z)\l'm') . (34) 

We also note some particular cases of expressions ()27|) and (jHUj) . First, consider 
the case s = Sq = z. This corresponds to the source and detector being oriented 
perpendicular to the surface of a slab. We then use Yim{z;z) = 6mo\/ (2/ + l)/47r to 
obtain 



^ ^ V(2/ + l)(2/^ + l) f d\ 
G(r,z;ro,z)= ^ / ——^exp[iq-{p-po)]{lO\K{q-,z-Zo)\lO). 

i,i'=o ^ J { ''^J 

(35) 

With the use of identity (IXl2|l given below, {lO\K{q; z)\l'0) can be expressed in terms 
of the associated Legendre functions of the first kind P[^{x) as 



V rni=—lm2=—l f M 

X (^rn#,)^^^Pl=^4^(V^,|/'m2)Pr MM)] ■ (36) 
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Next, consider the case q = 0. The operator k(0; z) describes one-dimensional 
propagation due to a planar source. We use d^^^,{0) = 6mm' to obtain 

// I fn \\T '\ Smm'[sgn{z)]^+^' exp(-|2;|/A^) 

{lm\K{0;z)\lm} = } {lm\il)^) —{i)^\lm). (37) 

3. 5. Plane wave and evanescent modes for the RTE 

Until now we considered solutions to the inhomogeneous RTE. However, the solution of 
boundary value problems requires knowledge of the general solution to the homogeneous 
equation. We seek such a solution in the form exp(— k ■ r) X]/m(^"^l'^)^»"(^i Upon 
substitution of this ansatz into the RTE with e = 0, we find that |k| must be the 
generalized eigenvalue (with the direction of k being arbitrary) and |c) the generalized 
eigenvector of the equation kR\c) = S|c), where the matrices R and S are defined by 
(|T!?jl - (fTHjl . Next we use (jHI) to express the spherical functions Yimis-jk) in terms of the 
spherical functions Yim{s;z) defined in a laboratory frame with the 2;-axis pointing in 
a selected direction. Then the general solutions to the homogeneous RTE (pj can be 
written as a superposition of the following modes: 

W(r,§) =exp j;il.(s;z)^^^PtpM4^(^.)(;|0^(M)) . (38) 

Here it is more convenient to use the notation for block eigenvectors |0„(M)) which 
were introduced in Section 13.2.11 The modes are labeled by the unit vector k 
(k ■ k = 1) whose polar angles in the laboratory frame are 6^^ and ip^^ and by the 
composite index /i = (M, n). We note that it is sufficient to consider only modes with 
positive eigenvalues Xmu (a^ ^ see Section I3.2.2|l due to the obvious symmetry 

^-k,-A/,n(-r,s) =/k,A/,n(r,s). 

However, the modes ()38|) with purely real vector k can not be used to construct a 
solution to a boundary value problem in a half-space or in a slab. This is because each 
mode is exponentially growing in the direction — k. Therefore, it is necessary to define 
evanescent modes with complex-valued vectors k. These modes are oscillatory in the 
X — y plane and exponentially decaying (or growing) in the z-direction. Namely, let 

k = -iXMnq. ± Z^MnQMniq) , (39) 

where q ■ z = 0. The polar angles of k are defined as follows: (p^ = (p^, cos(6'j^) = k ■ z = 
^^MnQMniQ) and sin(6'j^) = k ■ q = —iqXMn- Thus, we can write 9^^ = ir^qXMn), where 
the sine and cosine of the complex angle ?r(x) are given by ((221) ■ This gives rise to two 
kinds of evanescent modes: 



z I 



X dUiMqXMn)]{l\MM)) , (40) 
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X {-iy+"^di^_^[zT{qXMn)]{l\MM)) . (41) 

Here we have used the symmetry properties of rf^A/l^) under the transformation 
6^71 — 6 (which corresponds to change of sign of the factor cos(^)); hence the additional 
phase factor (—1)'+™ and the different sign of the index M in (j4H) . The phase factor 
(— 1)^^ is introduced for convenience. 

The plane-wave decomposition (j3(J|) can be equivalently rewritten as an expansion 
in terms of evanescent waves: 

(27r) 



G(r, s; ro, §o) = / 7^2^^^>^ §) ^-^(^0, -§o) , (42) 



where 

and the upper signs must be selected if z > zq while the lower signs are selected if z < zq. 
It can be easily verified that the expression (|42j) obeys the reciprocity condition. 

Now we make the following observation. Evanescent waves propagating in different 
directions can not, in principle, have identical angular dependence. In particular, by 
analyzing only the angular dependence of specific intensity in the plane z = (in 
infinite space), it is possible to tell if this intensity was produced by sources located 
to the left of the observation plane (in the region z < 0), or to the right. To 
demonstrate this point, we introduce vectors \ri^{q)) = IrjMnil)) with components 
{lm\r]Mn{q)) = d'-^]yj[iT{q\Mn)]{l\<Pn{M)) and a set of vectors obtained from \ri^{q)) 
by coordinate inversion: \fi^{q)) = V\ri^{q)). The evanescent modes fHU Il . pTI) can be 
written as 



/W(r,s) = e^^[ic^- p-Q^{q)z]Y^Yi^{s-i){lm\A{(\)\7]^{q)) , (44) 

Im 

/H(r,s) = e^^[ic^-p + Q,{q)z]Y,Yas-mrn\A\-mM)) , (45) 

Im 

{lm\A{q)\l'm') = 5w5mm' exp(-zmv9q)/ (46) 

Here ^(q) is a diagonal matrix. Note that ^(q) depends only on the direction of 
the real vector q, while \Ti^{q)) depends only on its length. In the case g = we 
have |?7^(0)) = and |r/^(0)) = Thus, the set {|r7^(0)), |?7^(0))} forms an 

orthonormal basis in the Hilbert space spanned by the eigenvectors of W, Ti. In the 
case g 7^ 0, the vectors {\ri^{q)), IVfiiq))} are no longer orthonormal. However, we 
believe on physical grounds that they still form a basis in Ti. Then there exists 



% We do not have a direct proof of this statement. However, in the opposite case, the boundary value 
problem would not be uniquely solvable. 
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ajual basis {K^ (<?)), |Cm(9))} such that {Cf,{q)\'nM)) = V' 1^71^(9)) = S^,^ and 

{CM)\vM)) = {Uq)\Uq)) = 0- 



Now assume that we have measured the specific intensity in the plane z = 0. 
Denote the two-dimensional Fourier transform of this function with respect to x and 
y by /o(q, s) = /im(q)^Zm(s; z). The expansion coefficients /im(q) are elements 
of a vector |/(q)). For every value of q, we can form a vector ^~^(q)|/(q)), since 
^(q) is invertible. If the sources are located to the left of the measurement plane, 
then ((^^(g)|A^^(q)|/(q)) = for every /i G In other words, the projection of 

^^^(q)|/(q)) onto the dual subspace spanned by |C/^(q')) is equal to zero. Analogously, 
if the sources are located to the right of the observation plane, {(^{q)\A~^ {q)\I (q)) = 0. 
Since the projection of a nonzero vector onto both subspaces can not be simultaneously 
zero, the angular dependence of the specific intensity in the plane z = carries 
information about the location of the source with respect to this plane. We emphasize 
that this analysis is only valid in the absence of boundaries. 

3.6. Boundary value problem 

Any solution to the RTE in a half-space or in a slab can be constructed as an expansion 
over modes P^ . p5|l with indices fi corresponding to only positive eigenvalues A^. This 
fact is crucial for application of the MRRF to the boundary value problem, which in 
radiative transport theory is formulated in the half-range of the angular variable. Now 
we demonstrate how it can be used to construct a solution to the boundary value problem 
posed on one or two planar interfaces. 

3.6.1. External source incident on a half-space 

Consider the RTE in the half-space z > 0. In this Section we assume that there are 
no internal sources in the medium, i.e., the RTE has a zero source term. The presence 
of external sources is expressed through an inhomogeneous boundary condition at the 
interface z = 



Here Iq{p,s) is the specific intensity evaluated at z = and Jinc(p, s) is the intensity 
incident from vacuum (the external source). The boundary condition (jTfj) is formulated 
in the half-range of the singular variable. 

The general solution to the RTE in the half-space z > can be written as a 
superposition of outgoing evanescent waves of the form ()44|) : 



where the unknown coefficients iq,/^ must be found from the boundary condition ()47|1 . 
Now we use expansion (|lSj) and expression to calculate Jo(p, s). Upon Fourier 



/inc(p,s) , if s ■ z > . 



(47) 




(48) 
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transformation of (|47p with respect to p, we arrive at the following equation: 

5^5^Vz™(s;z)(/m|^(q)|r7^(g))FW = /inc(q,§), if s ■ z > 0. (49) 

Im 

Next, we multiply both sides of Eq. HHlby Yif^,{s]z) and integrate over all directions 
such that s • z > 0. Note that integration in the right-hand side can be extended to 
all directions of s since /inc(ci5 s) is identically zero for s ■ z < 0. Thus, for a collimated 
narrow incident beam which crosses the boundary at p = po in the direction Sq, we 
obtain 

Y^^{lm\BA{q)Mq))F^+J = exp(-^q ■ po)Y;^{so; z) , (50) 
where matrix B is given by 



{lm\B\l'm')= [ Yil{s;z)Yi,^,{s;z)dh 

Pr{x)Pir{x)dx . (51) 



5mm' /(2/ + l)(2/' + l)(/-m)!(/'-m)! 



{I + m)\{l' + m)\ 







For a fixed value of q, fj50p is a set of linear equations of infinite size. In practice, this 
set must be truncated so that / < /max- Then the number of equations is 2N = (/max+1)^, 
where we have assumed for simplicity that /max is odd. But the number of unknowns 
F^J is only equal to N, since /i G Therefore, (j50p is formally overdetermined. 
However, not all equations in ()50|) are linearly independent. In fact, the rank of B is 
exactly equal to half of its size, which is a consequence of half-range integration in ()51|) . 
Therefore we come to the conclusion that (jKn|l is a well-determined system of equations 
with respect to the unknowns F^J . 

Numerically, the problem of solving ()5()j) can be solved in two different ways. A 
direct approach is to consider only equations in ()50|) which are linearly independent. 
This is achieved by only leaving equations in the system with / having the same parity 
as m, e.g., for a fixed m, / = \m\, \m\ + 2, \m\ + 4, . . ., with the restriction / < /max- 
Another approach is to seek the generalized Moore-Penrose pseudoinverse of (jSUj) . In 
this case eigenvectors of the truncated matrix B must be found numerically. If the size 
of B is even, half of its eigenvalues will be zero. Let be the eigenvectors of B with 
nonzero eigenvalues Oj/. Then the system of equations ()50|1 can be rewritten as 

a.Y^^mAmv,iq))FtJ = exp(-^q ■ Po) J2mim)Y;jso; z) . (52) 

Im 

We note that in the limit /max oo, the eigenvectors of B are known and are of simple 
form: (/m|^s) = Yi^{s; z) with the eigenvalues being unity for s-z > and zero otherwise, 
i.e., B is idempotent. 

The system can be simplified by the substitution Pj^J = f^f} exp(— zq-po). The 
coefficients f^f} are then independent of the source coordinate po. Another simplification 
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is achieved by noting that both A and B are diagonal in indices m and m'. In effect, the 
system fl^ must be solved once for each value of |q|; the dependence of the solution on 
the direction of q is trivial. If, in addition, the incident beam is normal to the interface 
(sq = z), the solutions do not depend on q at all. 

The additional computational complexity associated with solving the boundary 
value problem is then as follows. For every numerical value of the lengths of the vector 
q which is used in the expansion (jlH)) . a system of linear equations of size = (/max + 1)^ 
must be solved (the cost of diagonalization of B is negligibly small). Thus, consideration 
of boundary conditions adds significant computational complexity to the problem. This 
is a consequence of the fact that the rotation matrices exp[r(gA^) Jj,], unlike the matrix 
W, are not diagonal in m and m'. As a result, the system of equations ()52j) is not 
block diagonal and, in addition, g-dependent. However, the problem is easily solvable 
for /max < 100, which is, perhaps, more than is needed in any practical computation. 

3.6.2. External source incident on a slab 

The generalization of the mathematical apparatus developed in Section 13.6.11 to the 
case of the RTE in a finite slab is straightforward. Consider RTE in the slab < z < L. 
The external source is assumed to be incident from the left. Then the boundary 
conditions read 



/o(p, S) = /inc(P, s) , if s ■ z > , (53) 
/l(p,s) = 0, ifs-z<0, (54) 

where Jq and II are the specific intensities evaluated at the surfaces z = and z = L, 
respectively. The general solution inside the slab has the form 

Hr, s)=l ^Y^^ [^cSi^Hr, §) + i^i/ii(r, -§)] , (55) 

where F^^} and F^J are unknown coefficients. After some manipulations, we arrive at 
the following system of equations: 

{{lm\BAmnM))Ft^+^M-QM)L] {l,-m\BA\ci)\nM))F^:i] = 

exp(-iq- po)>^^(so;z) , (56) 
{exp [-Q,{q)L] {lm\BA{q)\fi,{q))FM + (/, -m\BA^ {ci)\^,{q))F^J} = . (57) 

This set of equations is the analog of ()5()|1 for the case of a finite slab. In the limit 
L oo one has F^J = and ()56|) coincides with ()50|). We note that for a fixed q, 
fl56|l . (jKTj) is a set of 2N linearly independent equations for 2N unknowns. The methods 
briefly discussed in Section 13.6.11 can be used to obtain the solution. 
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3. 6. 3. Internal source in a half-space 

Next, we consider an internal source in the half-space z > 0. Here we assume that 
there are no external sources. However, if this is not so, the solution can be obtained 
by simple superposition. 

Consider a point unidirectional source of the form e = 6{p — po)S{z — zo)6{s — Sq), 
where Zq > 0. The general solution in the region < z < Zq is written as 

/(r, s) = l s)/i+i(ro, -So) + F^Jl^^^^{r, s)] , (58) 

The second term in the square brackets in the right-hand side of the above expression 
can be interpreted as the surface term in the Kirchhoff-type formula for the Green's 
function (for the formulation of the Kirchhoff integral specific to the RTE see [23] or, 
for a more detailed derivation, [26]). The boundary condition at the interface z = 
reads: 

Jo(p,s) = 0, ifs-z>0. (59) 

The fact that the boundary condition is homogeneous reflects the fact that there are no 
external sources. The latter can be included by considering an inhomogeneous boundary 
condition of the type ()47|) . By analogy with Section rj.6.1l we immediately arrive at the 
following set of equations for the unknown coefficients Fq~^J: 

J^^^ [{lm\BA{q)\v,{q))F'tJ + V^A^^ -m\BA\q)Mq))I%{ro, -§0)] = . (60) 

Similarly to (|5(J|). this is a set of linearly independent equations with respect to 
unknowns. 



4. Numerics 



Now we illustrate the expressions obtained in Section for the RTE Green's function 
in an infinite medium with several numerical examples. We have computed the Green's 
function by truncating the series in (j^Hj) at /, /' < /^ax and using the reference frame 
in which the 2;-axis is aligned with the direction of the source, Sq. Correspondingly, 
the expression ()26p was used to compute the matrix elements of x- Note that in this 
expression the summation over M and j is finite; however, the summation over the 
modes is infinite and must be truncated. We have found empirically that, for each 
block B{M) of the matrix W, summation over = 500 eigenmodes (which corresponds 
to 1, 000 X 1, 000 matrices B(M)) is sufficient for all cases shown below. Further increase 
of does not change the result within double precision machine accuracy. The results 
start to deviate noticeably from those computed at A^ = 500 when A^ is taken to be 
smaller than ~ 100, especially when Zmax is relatively large. Further, we have used 
the Henyey-Greenstein model for the phase function, so that Ai = were < (7 < 1 
is a parameter. In all figures shown below we calculate the specific intensity /(r, s) 
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due to a point unidirectional source placed at the origin and illuminating in the z- 
direction. The distance from the source is measured in units of the transport free path, 
i* = l/[/ia + (1 — 9)f^s], which plays an important role in diffusion theory. 

It should be noted that numerical implementation of the formulas derived in this 
paper requires a degree of caution because the Green's function of the RTE is not 
square integrable with respect to both of its arguments, r and s. Therefore, one can not 
expect uniform point-wise convergence of the result with /max- Mathematically, this is 
manifested in the fact that the Bessel functions ki{x) that enter into ()25j) . (12611 diverge 
factorially for large orders: ki{x) oc /!! (/ ^ oo). This growth can not be compensated 
either by the Clebsch-Gordan coefficients, or by the eigenvector components (/m|?/'^) 
which decay, at best, exponentially (see discussion in Section I3.2.3|) . Therefore, 
p3|l . p5|l . ()26|l must be viewed as expressions defining the moments of the Green's 
function and the latter as a generalized function or a distribution. Nevertheless, in most 
practical situations, the spatial and angular dependencies of the Green's function can be 
approximated by smooth square-integrable functions by truncating at certain values 
of /max that provide desirable angular resolution. Computations are further facilitated 
by analytical subtraction of the ballistic component of the Green's function: 

Gfe(r,s;ro,So) = ^(s - So)^(R - So) ^^^^^f , R = r - Tq . (61) 
The corresponding ballistic contribution to Xb is 

{lm\xb{R; So)\l'm') = (5^^o V(2/ + 1)(2/- + i f'^^'j^f^ , (62) 

However, it is impossible to remove the singularities completely, and the remainder of 
such a subtraction still remains non-square integrable. 

The effect of subtraction of the ballistic term and convergence with /max for forward 
propagation is illustrated in Fig. |21 Here 6 is the angle between the direction of 
observation, s, and the positive direction of the 2;-axis: cos 6' = s • z. In the left 
column (plots a,c) we show the dependence of the specific intensity (with the ballistic 
term subtracted) on the maximum order of spherical functions /max- We assumed that 
convergence was reached when incrementing /max by 1 resulted in less than 0.1% relative 
change of the specific intensity in any direction. However, we emphasize again that this 
convergence is asymptotic. In the right column of images (b,d), we compare the angular 
dependence of the specific intensity for the maximum value of /max which was used in 
the graph to the left with and without the ballistic term. Note that the subtracted 
ballistic term can be added back analytically to the solutions obtained. In all figures 
shown below, the ballistic term is subtracted. 

Fig. ini illustrates the specific intensity for forward and backward propagation. 
Optical parameters were chosen to be close to those of typical biological tissues in the 
near infrared spectral region (see figure caption for details). The point of observation 
is placed at r = {0,0, z), where z is positive for forward propagation and negative for 
backward propagation, and 9 is defined in both cases as the angle between the vector 
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Figure 2. Angular dependence of the specific intensity for forward propagation at 
the distance z from the source. Left column (a,c): convergence with parameter Zmax- 
Right column (b,d): solid line shows the converged result obtained with subtraction 
of the ballistic term; dashed line: result obtained without subtraction of the ballistic 
term for the same /max- Top row (a,b): g = 0.5, Ha/l^s = 0.5, z = 20i*. Bottom row 
(c,d): g = 0.2, Ha/ij.^ = 0.01 and z = 10£*. 




Figure 3. Angular dependence of the specific intensity for forward (a) and backward 
(b) propagation obtained at /max = 21, 5 = 0.98 and Ha/l->'s = 6 • 10~^. The distance 
to the source z is assumed to be positive for forward propagation and negative for 
backward propagation. 
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s and the positive direction of the z-axis. It can be seen from the figure that the 
specific intensity in the backward direction is significantly smaller compared to that in 
the forward direction, even at relatively large source-detector separations {\z\ = 6i*). 
It can be also seen that the angular distribution of the specific intensity in the forward 
direction is more sharply peaked than that in the backward direction. This can be 
explained by noticing that backward propagation involves more scattering events than 
forward propagation of the same distance. 

Now we turn to the off-axis case. Here the source is still placed at the origin and 
illuminates in the positive z direction, while the point of observation is placed at a point 
r = {0,y,0). Below, we show two type of graphs. In the first case, the vector s is in 
the y — z plane, and its orientation is characterized by the angle a with respect to the 
positive direction of the 2;-axis. In the second case, s is in the x — y plane (perpendicular 
to So) and is characterized by the angle /3 with respect to the positive direction of the 
y-axis. The angles a and (3 (not to be confused with the Euler angles) are illustrated in 
Fig. 13 Note that a varies from to 2tt while f3 is restricted to the interval [0, vr] due to 
the obvious symmetry. 

In Fig. El we illustrate the specific intensity for highly forward-peaked scattering 
[g = 0.98) and the following three different ratios of fia/fJ's'- 6 ■ 10~^, 0.03 and 0.2. Note 
that the corresponding ratios of fJ^a/ f^'s^ where /i'^ = (1 — g)fis is the reduced scattering 
coefficient, are 0.003, 1.5 and 10, respectively. In the first case, the transport mean 
free path is mainly determined by scattering, while in the third case it is determined 
by absorption. The left column of images (a,c,e) illustrates the angular dependence of 
the specific intensity as a function of the angle a (vector s is in the y — z plane). The 
oscillations visible in Fig. El^e) is due to the non-square integrability discussed above. 
However, the values of the specific intensity at the region where the oscillations are 
visible are two to three orders of magnitude smaller than those at the peak. 

It is interesting to analyze the position of the maximum of the curves in Fig.|ni^a,c,e), 
ao- As the distance between the source and the detector increases, ao approaches 
tt/2. This corresponds to a vector s coinciding with the direction from the source 
to detector. However, for relatively small source-detector separations, ao is larger 




Figure 5. Angular distribution of specific intensity for off-axis propagation. 
Parameters: g = 0.98 and fJ.a/^^s = 6 • 10^^ (a,b), fJ-a/^J-s = 0.03 (c,d), ^Ma/lJ-s = 0.2 
(e,f). 



than 7r/2. The dependence of on the source-detector separation is illustrated in 
Fig. infa) for physiological parameters. The dependence of ao on the source-detector 
separation can be understood at the qualitative level. Indeed, at large separations, 
the angular distribution of the specific intensity is expected to be independent of the 
source orientation, with the maximum attained when s is aligned with the direction 
from the source to the detector. This corresponds to ao = vr/2. At smaller separations, 
the "photons" arrive at the detector locations along some "typical" (most probable) 
trajectories which are schematically illustrated in Fig. IHfb). We assume here that ao is 
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Figure 6. (a) Dependence of the position of maximum ao on the distance to 
the source, y, for physiological parameters: g — 0.98 and Ha/f^s = 6 • 10~^. (b) 
Schematic illustration of typical "photon trajectories" that correspond to maxima in 
graphs inta,c,e). 
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Figure 7. Convergence of the specific intensity with l^^x for g = 0.98, /ia//^s = 0.2 
and y = 22£*. 



determined by the angle at which the most probable trajectory crosses the y-axis. 

In Fig. |ni^b,d,f), the specific intensity is shown as a function of the angle f3 
(vector s is in the x — y plane). In this case, the maximum of the curves always 
corresponds to (3 = 0, which could be also inferred from the symmetry. We note that 



Lyi/3 = 0) = I, 



yz 



[a = vr/2), where the lower subscripts indicate the plane in which 
contains the vector s. 

The curves shown in Fig. El have a dynamic range of approximately 10^. A dynamic 
range of this magnitude was obtained due to the use of large values of /max- For smaller 
values of /max, the result can be grossly inaccurate and even negative. For example, 
in Fig. 13 we illustrate convergence with /max to one of the curves shown in Fig. El^e) . 
An accurate value of specific intensity at a ~ 7r/2 10~^ relative error) was obtained 
at /max = 39. Note that at /max = 10, the computed specific intensity is still grossly 
inaccurate. 
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5. Discussion 

The theoretical approach developed in this paper is, essentially, a spectral approach. 
Spectral methods have been studied extensively for the one-dimensional RTE [27]. 
However, in the 3D case these methods become very difficult to use. The substantially 
novel element of this paper is that we derive a usable spectral method for the full 
three-dimensional RTE with an arbitrary phase function and planar boundaries. The 
analytical part of the solution is of considerable complexity. However, this complexity 
is traded for the relative simplicity of the numerical part. In fact, we believe that we 
have reduced the numerical part of the computations to the absolute minimum which 
is allowed by the mathematical nature of the problem. 

This paper is limited to consideration of spatially-independent optical coefficients 
and phase functions. However, we note that the Green's function for a macroscopically 
homogeneous medium is of special interest, since it is used in linearized image 
reconstruction in optical tomography [28] and, more generally, in nonlinear image 
reconstruction based on the inversion of a functional series or the Newton-Kantorovich 
method [29]. Assuming the presence of only absorptive inhomogeneities in the medium, 
the linearized kernel of the integral equation of diffusion tomography has the form (in 
the slab imaging geometry) [28] 



r(pi, P2; r) = / G{pi, z = 0, Si = z; r, s)G(r, s; p2, z = L,S2 = z)d'^s , (63) 



where pi and p2 are the transverse coordinates of the source and detector, respectively, 
located on opposite surfaces of the slab, and G is the slab Green's function with constant 
absorption and scattering coefficients. One of the advantages of solutions obtained in 
this paper, compared to those based on discrete ordinates, is that the angular integral 
in the above formula can be evaluated analytically. 

This research was supported by the NIH under grant P41RR02305. 

Appendix A. Calculation of the integral (j24j) 

In this Appendix we evaluate the integral ()24|) for different choices of z. Written in 
component form, this integral reads 



We start with the case z = r. Then we have k ■ r = A;r cos%. We also notice that 
{lmi\ip ^) {ip oc 5mim2i so that the summation over mi and m2 can be replaced by 




{lmi\%lj^,){%ljy\l'm2) 
1 + ik\n 
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summation over a single index M which runs from — f to T, where I = min(/,/'). Then 
flA.l|) can be rewritten as 

where / is the angular part of the integral (the list of formal arguments of I is omitted): 
1 = j ^^^^^'^^^'^'f'^ exp[^(m' - m)ip^] exp{tkr cos^k)4Af(^k)^!n'A/(^k) • (A.3) 

The integral over y^j^ is evaluated immediately with the result 2Tr6mm'- Integration over 
9^^ requires expanding the exponent in the integrand as 

oo 

exp(zA;r cos^^) = $^^''(2^ + l)jL(^^)4o(^k) ' (A-4) 

L=Q 

where JlIx) are the spherical Bessel functions of the first kind, and using the following 
formula (see Ref. [25], Sect. 4.11.2, formula 8, and the symmetry properties of d-functions 
given in Sect. 4.11, formula 1 of the same reference): 

o/ 1 \m~M 

where C^jlmlj2m2 Clebsch-Gordan coefficients. Taking into account that C/"^;, _^ 

is nonzero only for |/ — /'| < L < / + we obtain 

l+V 

L=\l-l'\ 

Next, we substitute this result into ()A.5j) and, after some rearrangement, arrive at 

{lm\x{r;m'rn') = '^-;^^"^ ^ (-1)" E .^CS,.,_™C^Sa_m 

^/c^/C^// '—^ 

^ M=-l L=\l-V\ 



To evaluate the radial integral, we exploit the symmetry properties of the above 
expression. First, we notice that Cf'^i, = (— l)'+''+'^C/'l°^,f j^^, while {lM\'(p^) does 
not depend on the sign of M. Thus, the addition of terms with positive and negative 
values of M in the above formula (for M 7^ 0) gives zero unless l+l'+L is even. Likewise, 
in the case M = 0, C/"q°, q = unless the above sum of indices is even. Correspondingly, 
the only nonzero contributions to the sum over L corresponds to L = \l — 1'\ +2j , where 
the index j runs from to /. Next, we use the symmetry property of the eigenvectors 
discussed in Section 13.2.21 This property allows one to limit summation over the 
eigenvector indices fi to only the values corresponding to positive eigenvalues while 
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simultaneously replacing the factor 1/(1 + ik\^) by 1/(1 + ik\^) + (— 1)^+^'/ (1 — ik\^). 
Thus, we obtain 



Y^{lM\ij,){ij,\l'M) J, (A. 



X 

where J is the radial integral given by 

^"io (2^^i'-"i+^^-^^"^ ttpai ■ ^^-^^ 

The parity of the Bessel functions in the above integral is the same as that of / + 
Therefore, the integrand is an even function of k for all values of the indices, and the 
integral can be extended to — oo and calculated by residues. The result is 

J = vrr(i'-''i+^^-)A;3fc|,_,|+2,(r/A^) . (A.IO) 

Upon substitution of this result into ()A.8|1 . we obtain the formula (j^ . 

In the case z = Sq the dot product k ■ r can not be written as kr cos6^. Therefore, 
the exponent in the angular integral I is expanded as 

exp(zk • r) = 47rJ2^''JL{kr)YLM'{k So)Ylj,,,{r; Sq) . (A. 11) 

LM' 

We further take advantage of the identity 

YLM'i0i,,Vi,) = (-l)"^V4vr/(2L + l)4/'(^k) exp(zMVk) (A.12) 



to transform the angular integration to the general form ()A.6|) . Note that azimuthal 
integration results in a factor of 6M'm and thus removes summation over M'. The final 
result for / is 

^ = ^^^= E ^"jl(^0>^£™(?;so)<.J:,,o^,:;C.,o, (A.13) 

+ J- L=\l-l'\ 

where we have also used C^%^,_m = (-1)'"*V(2^ + l)/(2^' + 1)Cwl,o- The radial 
integration, and the symmetry considerations explained above, remain without change. 
Substitution of ()A.13|) into ()A.2|) and subsequent radial integration leads to the formula 

[E 
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Appendix B. Calculation of the integral (|28|) 

Integral ()28|1 . written in components, reads 

(H^(q;^)|/W) = --P[-(--"^)] f y-mm^AMM (b.i) 

Here I is tlie integral over /c^: 

where we have used the notations introduced in Section IH.2.11 for block eigenvectors 
|0n(M)). The angle 6 is defined by fl^ in Section The Wigner d-functions can be 
written in terms of cos 9 as 

m-M| |7Ti + M| 

where ^rnM = 1 if m < M and ^ = (-1)™+^^ if m > M, 



1^- 


m - M\/2 - 


m + M\/2)\{l + 


m + M\/2- 


m + M\/2)\ 


{1 + 


m - M\/2 - 


m + M\/2)\ {I - 


m + M\/2- 


m + M\/2)\ 



and Pi"'^^(x) in expression ()B.3j) are Jacobi polynomials with s = I — \m — M\/2 — \m + 
M\/2, u=\m- M\ and = |m + M|. 

The integrand in ()B.2|) is not, in general, an analytic function of k^- However, the 
expression for the Green's function contains a summation over M. It can be shown 
explicitly that the combination 

dlMmlMiO) + dl^,,{e)dl_,M (B.5) 

contains only even powers of the factor ^Jk1 + if / + /' is even and only odd powers of 
the same factor if / + /' is odd (a general proof of this statement is available but omitted). 
Taking into account the factor \J + k1[l — (— l)'"*"''] in the right-hand side of ()B.2|1 . 
we arrive at the conclusion that the integrand becomes analytic after addition of terms 
with positive and negative values of M. Note that the eigenvectors and eigenvalues do 
not depend on the sign of M and the above consideration applies to the case M = 0. 
Consequently, one can evaluate ()B.2|) by residues choosing a branch of the complex- 



valued function \/k1 + arbitrarily. 

The integrand of ()B.2|) has simple poles at kz = ±i^/q'^ + Taking account 

of these poles leads to the following expression: 



[sgn(z)] 



l+l'+m+m' 



Xmu exp 



-a/1 + {q>^Mny\z\/XMn 



'1 + {qXMny 
X dLM{^^i<l^Mn)]di,'M[i'r{q\Mn)] 



(B.6) 
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Substitution of ()B.6|) into ()B.1|) leads to an expression which is equivalent to (j2Hl)- 

We note that the integrand of ()B.2|) has another set of poles. Namely, these are 
poles of the functions d[^j^[9{kz)] at kz = ±iq- These poles are of a purely geometrical 
nature. We have calculated analytically the contributions of these poles to the Green's 
function to the few lowest orders in I, I', and found that they cancel each other. However, 
we do not have a general proof of such cancellation to all orders. On the other hand, it is 
clear that if these poles could contribute to the plane- wave decomposition of the Green's 
function, the result would not satisfy the RTE since the matrix W is bounded and has 
no infinite eigenvalues. To confirm the validity of the obtained analytical expression, 
we have computed / numerically by the fourth-order Simpson rule for a model set of 
parameters. Then we used this result to compute the Green's function for the particular 
case s = So = z. The result coincided with the one predicted by formula p5|) with 
machine accuracy in double precision. 
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